Inferência Causal - 2024/2
A inferência causal é uma área de estudo de interesse em várias áreas do conhecimento, como por exemplo a saúde, a economia.
No entanto, a complexidade dos problemas a serem resolvidos e o volume de dados a ser usado é cada vez maior.
Na inferência tradicional, podemos nos utilizar de técnicas de Machine Learning para contornar esse problema, uma vez que estes modelos são bons para lidar com problemas complexos e não exigem premissas fortes, como por exemplo linearidade.
Como podemos então utilizar Machine Learning para medir efeito causal?
É isto que propõe a técnica de Double Machine Learning (DML)!
Suponha que queremos medir o efeito do preço de nosso sorvete nas nossas vendas, controlando pelos confundidores: temperatura, custo e dia da semana.
A princípio suponha que usemos um modelo linear para remover o viés de nosso efeito causal:
\[ Vendas_i=\alpha+\tau preço_i+\beta_1 temp_i + \beta_2 custo_i +\beta_3 dsem_i+\epsilon_i \]
Só estamos interessados no parâmetro \(\tau\), mas como a estimação dos parâmetros é feita simultaneamente, precisamos nos preocupar na estimação dos demais parâmetros, pois afetarão a estimação do efeito causal.
Por exemplo, a relação de \(temp\) provavelmente não é linear, pois em dias muito frios as pessoas não costumam sair, o que afetará nossas vendas, mas o mesmo pode ser dito para dias muito quentes, em que está tão quente que as pessoas preferem ficar em casa, no ar condicionado.
Dado isso, um modelo mais preciso seria o seguinte:
\[ Vendas_i=\alpha+\tau preço_i+\beta_1 temp_i^2 + \beta_2 custo_i +\beta_3 dsem_i+\epsilon_i \]
Ou seja, note que quanto maior o número de variáveis mais complexo é definir essas interações. Devemos colocar quais polinômios, quais interações?
É um processo muito complexo que exige um estudo aprofundado de conhecimento de domínio.
Considere o seguinte modelo inicial,
\[ Y_i = \tau T_i + \beta X_{i} + e_i \]
Vamos ajustar duas regressões lineares:
\[ Y_i \sim X_{i} \\ T_i \sim X_{i} \]
A partir desses modelos obtemos dois resíduos, que chamaremos de \(r_Y\) e \(r_T\), respectivamente.
Agora, considere a seguinte regressão:
\[ r_Y \sim r_T \]
O teorema de Frisch-Waugh-Lovell nos garante que esta regressão retorna o parâmetro \(\tau\)!
Call:
lm(formula = resid(fit_y) ~ resid(fit_t))
Residuals:
Min 1Q Median 3Q Max
-3.7644 -0.6691 -0.0030 0.6658 3.7492
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.102e-16 1.005e-02 0.00 1
resid(fit_t) 1.997e+00 2.014e-02 99.17 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.005 on 9998 degrees of freedom
Multiple R-squared: 0.4959, Adjusted R-squared: 0.4958
F-statistic: 9834 on 1 and 9998 DF, p-value: < 2.2e-16
Que maravilha! Então podemos substituir as regressões intermdiárias por modelos não-paramétricos e agora não temos mais que nos preocupar com os demais parâmetros.
Mas peraí, podemos fazer isso?
Note que o teorema de Frisch-Waugh-Lovell (FWL) vale apenas para o cenário cujas relações são lineares.
Para nossa felicidade pessoas incríveis já extenderam esses resultados.
A princípio o resultado do FWL foi extendido para métodos “tradicionais” não-paramétricos, como Lasso (Bickel et al, 1998).
Os resultados teóricos nesse momento se limitavam a estimadores que tomassem valor num conjunto de Donsker, o que eliminava muitos dos novos métodos.
Mais recentemente, conseguiram expandir os resultados para métodos mais novos como florestas, gradient boosting e outros (Chernozhukov et al. - 2014).
Note que no teorema de FWL temos duas etapas intermediárias, que buscam estimar essas duas informações:
\[ E[Y|X] \\ E[T|X] \]
O que o DML propõe é utilizarmos métodos não paramétricos para estimá-las.
Com o uso de modelos não paramétricos teríamos uma maneira de contornar as restrições de ter que nos atentar em todos parâmetros, deixando isso para esses modelos intermediários e focar nossa atenção em estimar o parâmetro causal.
Com isso, bastaríamos fazer a regressão linear dos resíduos desses modelos e teríamos o efeito causal \(\tau\).
Vale ressaltar que quando analisado o efeito causal, não estamos falando de um efeito linear de fato, mas sim de um efeito linear local, sendo necessario cuidado ao analisá-lo para nao generalizá-lo para todo o intervalo real.
Como temos modelos de machine learning envolvidos, temos que nos preocupar com a questão do overtfitting.
No caso do debiased/orthogonal machine learning, temos dois modelos de regressão (\(Y\sim X\) e \(T\sim X\)) com os quais precisamos nos preocupar.
Falemos um pouco sobre as consequências de sobreajustar \(Y\sim X\) e \(T\sim X\). Lembremos que o objetivo do DML é ajustar um modelo com os resíduos dos modelos citados.
Nesse caso, os resíduos \(r_Y\) serão menores do que deveriam.
Ou seja, além de desconsiderar a relação entre \(Y\) e \(X\), uma parte da relação entre \(Y\) e \(T\) será perdida. Consequentemente, os resíduos estarão enviesados em torno do zero, já que haverá pouco a representar para a regressão final que utilizará os resíduos.
Se esse modelo for sobreajustado, então \(T\) será explicado mais do que deveria por \(X\), ou seja, a variância de \(T\) será explicada demais por \(X\).
Considerando o modelo final com o resíduo \(r_T\), não haverá muita variância restante nesses resíduos e então a variância do estimador final será muito alta, já que será como se todos os valores de \(r_T\) fossem iguais.
Como possivelmente solucionar esses problemas?
Utilizaremos cross-validation!
Ao invés de utilizar a mesma amostra para estimar e obter os resíduos, dividiremos os dados em k partes iguais.
O modelo será treinado com k-1 amostras e os resíduos serão calculados utilizando a k-ésima parte, esse processo será realizado para todas as amostras.
Desse modo, mesmo com overfitting, os resíduos de cada uma das regressões não serão enviesados, já que usamos uma amostra diferente da qual o modelo foi treinado.
Por fim, a estimação final é obtida combinando as predições de todas as k-1 amostras.
Para exemplificar o DML em ação, simularemos dados.
gen_data <- function(n, d, p, delta, base) {
X <- matrix(rnorm(n * d), nrow = n, ncol = d)
D <- rbinom(n, 1, p)
y0 <- base - X[, 1] + rnorm(n, mean = 0, sd = 0.1)
y1 <- delta + base - X[, 1] + rnorm(n, mean = 0, sd = 0.1)
y <- y1 * D + y0 * (1 - D)
return(tibble::tibble(y = y, D = D, X = X))
}
n <- 3000 # n samples
d <- 10 # n features
delta <- 7.0 # treatment effect
base <- 0.3 # baseline outcome
# Set the random seed for reproducibility
set.seed(125)
# Generate RCT data
modelData <- gen_data(n, d, 0.2, delta, base) |>
dplyr::mutate(DF = factor(D, levels = c(0, 1)))Agora faremos os dois modelos intermediários:
train_control <- caret::trainControl(
method = "adaptive_cv",
number = 10,
search = "random",
verboseIter = FALSE
)
rfResponseModel <- caret::train(
y ~ . - D - DF,
method = "ranger",
tuneLength = 10,
data = modelData,
verbose = FALSE,
trControl = train_control
)
rfTreatmentModel <- caret::train(
DF ~ . - y - D,
method = "ranger",
tuneLength = 10,
data = modelData,
verbose = F,
trControl = train_control
)Agora vamos calcular os resíduos de ambos os modelos:
# Predict the response in dataset 1 (2) using model 2 (1).
responsePredictions <- predict(rfResponseModel,modelData)
# Do the same for the treatment model
treatmentPredictions <- as.numeric(predict(rfTreatmentModel,modelData))-1
# Calculate the treatment residuals
treatmentResiduals <- modelData$D - treatmentPredictions
# Calculate the response residuals
responseResiduals <- modelData$y - responsePredictionsPor fim, faremos o modelo final com os resíduos:
Call:
lm(formula = responseResiduals ~ treatmentResiduals)
Residuals:
Min 1Q Median 3Q Max
-2.2137 -0.2548 -0.0854 0.0973 7.6734
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.15045 0.01604 -71.72 <2e-16 ***
treatmentResiduals 6.18765 0.03700 167.25 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7918 on 2998 degrees of freedom
Multiple R-squared: 0.9032, Adjusted R-squared: 0.9032
F-statistic: 2.797e+04 on 1 and 2998 DF, p-value: < 2.2e-16